#ifndef LBWRAPPER_HPP
#define LBWRAPPER_HPP

#include "LBWrapper.h"
#include <map>
#include <string>
#include <sstream>

////////////////////////////////////////// LBMBase //////////////////////////////////////////
template<typename T, template<typename U> class Lattice>
LBMBase<T, Lattice>::~LBMBase()
{  
	if(lattice) delete lattice;  
	if(boundaryCondition) delete boundaryCondition; 		
}

template<typename T, template<typename U> class Lattice>
void LBMBase<T, Lattice>::InitialLattice(T nx, T ny, T nz)
{
#ifdef PARALLEL_MODE_MPI
		this->lattice = new MultiBlockLattice3D<T, Lattice>(createRegularDataDistribution(nx, ny, nz));
#else
		this->lattice = new BlockLattice3D<T, Lattice>(nx, ny, nz);
#endif // PARALLEL_MODE_MPI
}

template<typename T, template<typename U> class Lattice>
void LBMBase<T, Lattice>::InitialBC()
{
	this->boundaryCondition = createLocalBoundaryCondition3D(*(this->lattice));
	//this->blockParamoundaryCondition = createInamuroBoundaryCondition3D(*(this->lattice));
	//this->boundaryCondition = createInterpBoundaryCondition3D(*(this->lattice));
}

template<typename T, template<typename U> class Lattice>
void LBMBase<T, Lattice>::ApplyNeumannBoundaryConditions()
{
	std::cout << "Apply Neumann Boundary" << std::endl;
	//LBM Global Data
	int nx = this->lattice->getNx();
	int ny = this->lattice->getNy();
	int nz = this->lattice->getNz();
	int Q = 15;

	//Parallel local data
	int iBlock = singleton::mpi().getRank();
	BlockParameters3D const &blockParam = lattice->getMultiData().getBlockParameters(iBlock);
	BlockCoordinates3D const &envelope = blockParam.getEnvelope();
	BlockLattice3D<T,Lattice> *localLattice = lattice->getBlockLattices()[iBlock];

	int nxLocal = blockParam.toLocalX(nx);
	int nyLocal = blockParam.toLocalX(ny);
	int nzLocal = blockParam.toLocalX(nz);

	int BoundaryPlate2nLocal = blockParam.toLocalZ(0);
	int BoundaryPlate2pLocal = blockParam.toLocalZ(nz-1);

	for (int iX=0; iX != nx; ++iX)
	{
		for (int iY=0; iY != ny; ++iY)
		{
			Cell<T, Lattice>& bmCell = lattice->get(iX, iY, 0);
			Cell<T, Lattice>& cellP = lattice->get(iX, iY, 0+1);
			Cell<T, Lattice>& cellPP = lattice->get(iX, iY, 0+2);

			T vel[3] = {T(), T(), T()};

			bmCell.computeU(vel);
			vel[0] = vel[1] =  0.;

			T rho = (4*cellP.computeRho() - cellPP.computeRho()) / 3.;
			//T rho = cellP.computeRho();

			bmCell.defineRho(rho);

			Cell<T, Lattice>& bpCell = lattice->get(iX, iY, nz-6);
			Cell<T, Lattice>& cellM = lattice->get(iX, iY, nz-6 -1);
			Cell<T, Lattice>& cellMM = lattice->get(iX, iY, nz-6 -2);


			vel[0] = T(); vel[1] = T(); vel[2] = T();

			bpCell.computeU(vel);
			vel[0] = vel[1] =  0.;

			rho = (4*cellM.computeRho() - cellMM.computeRho()) / 3.;
			//rho = cellM.computeRho();

			bpCell.defineRho(rho);

			for(int iGhost = 1; iGhost <= 4; ++iGhost)
				lattice->get(iX, iY, nz-6 + iGhost).defineRho(rho);

		}
	}

	// Z- direction von Neumann BC
//	if( (BoundaryPlate2nLocal > envelope.z0 && BoundaryPlate2nLocal < envelope.z1))
//	{
//		for(int iX = envelope.x0; iX != envelope.x1 + 1; ++iX)
//		{
//			for(int iY = envelope.y0; iY != envelope.y1 + 1; ++iY)
//			{
//				Cell<T, Lattice>& bCell = localLattice->get(iX, iY, BoundaryPlate2nLocal);
//				Cell<T, Lattice>& cellP = localLattice->get(iX, iY, BoundaryPlate2nLocal+1);
//				Cell<T, Lattice>& cellPP = localLattice->get(iX, iY, BoundaryPlate2nLocal+2);
//
//				bCell.defineRho( (4*cellP.computeRho() - cellPP.computeRho())/3. );
//			}
//		}
//
//	}
//
//	// Z+ direction von Neumann BC
//	if( (BoundaryPlate2pLocal > envelope.z0 && BoundaryPlate2pLocal < envelope.z1))
//	{
//		for(int iX = envelope.x0; iX != envelope.x1 + 1; ++iX)
//		{
//			for(int iY = envelope.y0; iY != envelope.y1 + 1; ++iY)
//			{
//				Cell<T, Lattice>& bCell = localLattice->get(iX, iY, BoundaryPlate2pLocal);
//				Cell<T, Lattice>& cellM = localLattice->get(iX, iY, BoundaryPlate2pLocal-1);
//				Cell<T, Lattice>& cellMM = localLattice->get(iX, iY, BoundaryPlate2pLocal-2);
//
//				bCell.defineRho( (4*cellM.computeRho() - cellMM.computeRho())/3. );
//			}
//		}
//
//	}

}


template<typename T, template<typename U> class Lattice>
void LBMBase<T, Lattice>::ApplySphereCynliderBC(T cRadius, T drift)
{
	std::cout << "Apply Cylinder" << std::endl;

	int iBlock = singleton::mpi().getRank();
	BlockParameters3D const& blockParam = lattice->getMultiData().getBlockParameters(iBlock);
	BlockCoordinates3D const& bulk = blockParam.getBulk();
	BlockLattice3D<T,Lattice> *localLattice = lattice->getBlockLattices()[iBlock];

	int centerX = lattice->getNx() / 2;
	int centerY = lattice->getNy() / 2;
	int centerZ = lattice->getNz() / 2;

	std::cout << centerX << " " << centerY << " " << centerZ << std::endl;

//	if(rank == 0)
//	{
//		std::ofstream debug("debug_0.txt");
//		debug << "[ " << rank << " ]-" << centerX << " " << centerY << " " << centerZ << " " << std::endl;
//		debug << "[ " << rank << " ]-" << bulk.x0 << " " << bulk.x1 << std::endl;
//		debug << "[ " << rank << " ]-" << bulk.y0 << " " << bulk.y1 << std::endl;
//		debug << "[ " << rank << " ]-" << bulk.z0 << " " << bulk.z1 << std::endl;
//	}
//	if(rank == 1)
//	{
//		std::ofstream debug("debug_1.txt");
//		debug << "[ " << rank << " ]-" << centerX << " " << centerY << " " << centerZ << " " << std::endl;
//		debug << "[ " << rank << " ]-" << bulk.x0 << " " << bulk.x1 << std::endl;
//		debug << "[ " << rank << " ]-" << bulk.y0 << " " << bulk.y1 << std::endl;
//		debug << "[ " << rank << " ]-" << bulk.z0 << " " << bulk.z1 << std::endl;
//	}

	//std::map< int, std::map<int, T> > phi2D;

	//std::vector< std::vector< std::vector<T> > > phi2D(51, std::vector< std::vector<T> >(51, std::vector<T>(201, 0)));

	for(int iX = bulk.x0; iX <= bulk.x1; ++iX)
	{
		for(int iY = bulk.y0; iY <= bulk.y1; ++iY)
		{
			for(int iZ = bulk.z0; iZ <= bulk.z1; ++iZ)
			{
				sMathLib::Vector3 rVector(iX-centerX, iY-centerY, 0.);

				T r(rVector.abs());

				T dist = r - cRadius;

				T phi;

				if( dist < -drift )
					phi = 0;
				else if( dist > -drift && dist < 0)
					phi = 1.-(dist + drift)/(drift + drift);
				else
					phi = 1.;

					//phi2D[iX][iY][iZ] = phi;

				int local_X = blockParam.toLocalX(iX);
				int local_Y = blockParam.toLocalY(iY);
				int local_Z = blockParam.toLocalZ(iZ);

				Cell<T,Lattice>& cell = localLattice->get(local_X, local_Y, local_Z);


				const T uB[3] = {0., 0., 0.};

				const T uSqr = 0.;

				const T density = cell.computeRho();

				for (int iPop=0; iPop != Lattice<T>::q; ++iPop)
				{
					cell[iPop] = (1. - phi) * cell[iPop] + phi * cell.computeEquilibrium(iPop, density, uB, uSqr);
				}

			}
		}
	}

//	std::ofstream phi("phi.dat");
//	phi << "variables = x, z, phi" << std::endl;
//	phi << "zone i = " << bulk.x1 - bulk.x0 + 1 << ", j = " << bulk.y1 - bulk.y0 + 1 << " , f = point" << std::endl;
//
//	for(int iY = bulk.y0; iY <= bulk.y1; ++iY)
//	{
//		for(int iX = bulk.x0; iX <= bulk.x1; ++iX)
//		{
//			phi << iX << " " << iY << " " << phi2D[iX][iY][10] << std::endl;
//		}
//	}
//
//	phi.close();



}

////////////////////////////////////////// LBMBase //////////////////////////////////////////

////////////////////////////////////////// IBMLin //////////////////////////////////////////
IBMLin::IBMLin()
{
#ifdef PARALLEL_MODE_MPI	
#endif
}

IBMLin::~IBMLin()
{
	if(bulkDynamics) 
		delete bulkDynamics;
}

void IBMLin::InitialDynamics_BC_Geom(T omega, sMathLib::VectorX<int, 3> domain, BoundaryInput<T>& bInput)
{
	this->bulkDynamics = new BGKdynamics<T, D3Q15Descriptor>( omega, instances::getBulkMomenta<T, D3Q15Descriptor>() );
		
//////////////////////////// Initial BC ////////////////////////////
	InitialBC();
//////////////////////////// Initial BC ////////////////////////////

	initGeometry<T, D3Q15Descriptor>( *(this->lattice), domain, omega, *(this->bulkDynamics), *(this->boundaryCondition), bInput);
	//continueGeometry<T, D3Q15Descriptor>( *(this->lattice), domain, omega, *(this->bulkDynamics), *(this->boundaryCondition), bInput);
}

void IBMLin::PrepareToCalculateInteraction(double externalForce[3]) 
{
	this->lattice->forAll(InitExternalForce<T, D3Q15Descriptor>(externalForce));
}
	
void IBMLin::PrepareToCalculateInteraction()
{
}
	
void IBMLin::ClearAllSphere()
{ 
	this->intCalu.clearAll();
}
	
void IBMLin::AddSphereImpl(const NJRvector3d& pos, const NJRvector3d& vel, const NJRvector3d& aVel, double r, double rotateMatrix[3][3])
{
	this->intCalu.addDO(pos, vel, aVel, r);
}

void IBMLin::CalculateSphereFluidInteraction(T LatticeNu, int iT) 
{		
#ifdef PARALLEL_MODE_MPI
	intCalu.parallelCalculateInteraction(*lattice, LatticeNu, iT);
#else
	intCalu.calculateInteraction(*lattice, LatticeNu, iT);
#endif // PARALLEL_MODE_MPI
}

void IBMLin::PrepareToNextStep()
{
}

T IBMLin::getVolumeFraction(int iX, int iY, int iZ, int np)
{
	return T();
}

void IBMLin::recordFluidEffect(int iT, int recordStep, std::vector<T>& scalingFactor)
{
	this->intCalu.recordFluidEffect(iT, recordStep, scalingFactor);
}

////////////////////////////////////////// IBMLin //////////////////////////////////////////

////////////////////////////////////////// LBWrapper //////////////////////////////////////////
template <class IBMPolicy>
LBWrapper<IBMPolicy>::LBWrapper(int argc, char* argv[], sMathLib::VectorX<int, 3> N_, double tau, const string& dir, T vtkLogStep, BoundaryInput<T>& bInput)
:_tau(tau)
{
///////////// OLB initialize /////////////
	olbInit(&argc, &argv, true);
  singleton::directories().setOutputDir(dir);
///////////// OLB initialize /////////////
  
	this->_domain = N_;
	this->_omega = 1. / _tau;
	
	this->iT = 0;
	this->_vtkSave = vtkLogStep;
	
	this->pbc = 0;
	this->LBM_periodic = false;	


//	std::ostringstream ostr;
//	ostr << argv[6];
//	std::istringstream

///////////// Initialize Lattice /////////////
	IBMPolicy::InitialLattice(_domain.x(), _domain.y(), _domain.z());
	IBMPolicy::InitialDynamics_BC_Geom(this->_omega, this->_domain, bInput);
///////////// Initialize Lattice /////////////
	
	std::cout << IBMPolicy::lattice->getNz() << endl;
	
	IBMPolicy::intCalu.setTau(tau);

	IBMPolicy::intCalu.createVolumeFractionField(*(IBMLin::lattice));

//////// Test member
	this->_pressure = bInput.getBoundaryValue("Z-Left-hand-side");
	IBMLin::intCalu.assignComputationalDomain(*(IBMLin::lattice));

	IBMLin::intCalu.createVolumeFractionField(*(IBMLin::lattice));
}

template <class IBMPolicy>
LBWrapper<IBMPolicy>::~LBWrapper() 
{
	
}

template <class IBMPolicy>
void LBWrapper<IBMPolicy>::SetPeriodicBoundaryConditions(const Boundary* pbc) 
{ 
	this->pbc = pbc; 
	this->LBM_periodic = pbc->GetSwitch(0) || pbc->GetSwitch(1) || pbc->GetSwitch(2);
}

template <class IBMPolicy>
void LBWrapper<IBMPolicy>::AddSphere(const NJRvector3d& pos, const NJRvector3d& vel, const NJRvector3d& aVel, double r, double rotateMatrix[3][3])
{
		IBMPolicy::AddSphereImpl(pos, vel, aVel, r, rotateMatrix);
}

//template <>
//void LBWrapper<IBM2008>::AddSphere(const NJRvector3d& pos, const NJRvector3d& vel, double r, double rotateMatrix[3][3])
//{
//
//}

template <class IBMPolicy>
int LBWrapper<IBMPolicy>::GetNumberOfSphere() const 
{ 
	return IBMPolicy::intCalu.GetNumberOfDO(); 
}

template <class IBMPolicy>
const NJRvector3d& LBWrapper<IBMPolicy>::GetSphereForces(int i)
{
	return IBMPolicy::intCalu.getForces(i);
}

template <class IBMPolicy>
const NJRvector3d& LBWrapper<IBMPolicy>::GetSphereTorque(int i)
{
	return IBMPolicy::intCalu.getTorque(i);
}
	
template <class IBMPolicy>
void LBWrapper<IBMPolicy>::ModifyBCVelocity(T base_vel)
{
}

template <class IBMPolicy>
void LBWrapper<IBMPolicy>::NextStep()
{	
	IBMPolicy::lattice->collideAndStream(this->LBM_periodic);		
}

template <class IBMPolicy>
void LBWrapper<IBMPolicy>::NextStep(int iT, ibm::LBMUnit<int, T>* cu)
{
	IBMPolicy::lattice->collideAndStream(this->LBM_periodic);

	
	if(iT % (int)_vtkSave == 0)
	{
		DrawField(iT, 0, cu);
	}


}

template <class IBMPolicy>
void LBWrapper<IBMPolicy>::DrawField(int iT, int interval, ibm::LBMUnit<int, T>* cu)
{
	WriteVTK(*(IBMPolicy::lattice), iT, cu);
}

template <class IBMPolicy>
void LBWrapper<IBMPolicy>::AssignBottomPressure(T pre)
{
	applyBoundaryPressure<T, D3Q15Descriptor>( *(IBMPolicy::lattice), pre);
}

////////////////////////////////////////// LBWrapper //////////////////////////////////////////

#endif //LBWRAPPER_HPP
